// Copyright 2012 Sergio, Raul, Dorival. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

// STL
#include <iostream>
#include <cmath>

// Mechsys
#include <mechsys/dem/domain.h>
#include <mechsys/mesh/structured.h>

using std::cout;
using std::endl;

int main(int argc, char **argv) try
{

    //////////////////////////////////////// 3D: structured //////////////////////////////////////////////////
    
    {
        /*
                          4----------------7
                        ,'|              ,'|
                      ,'  |            ,'  |
                    ,'    | -60   -10,'    |
                  ,'      |        ,'      |
                5'===============6'        |
                |         |      |    -40  |
                |    -30  |      |         |
                |         0- - - | -  - - -3
                |       ,'       |       ,'
                |     ,' -20     |     ,'
                |   ,'        -50|   ,'
                | ,'             | ,'
                1----------------2'
        */
        double Lx=20.0, Ly=10.,  Lz=10.;
        double          ly=0.5,  lz=0.5;
        Array<Mesh::Block> blks(9);

        // lower
        blks[0].Set (/*NDim*/3, /*Tag*/-1, /*NVert*/8,
                     -1.,  0.0, 0.0, 0.0,
                     -2.,   Lx, 0.0, 0.0, 
                     -3.,   Lx,  Ly, 0.0, 
                     -4.,  0.0,  Ly, 0.0,
                     -5.,  0.0, 0.0,  Lz,
                     -6.,   Lx, 0.0,  Lz, 
                     -7.,   Lx,  Ly,  Lz, 
                     -8.,  0.0,  Ly,  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[1].Set (/*NDim*/3, /*Tag*/-2, /*NVert*/8,
                     -1.,  0.0, Ly + 0.0, 0.0,
                     -2.,   Lx, Ly + 0.0, 0.0, 
                     -3.,   Lx, Ly +  ly, 0.0, 
                     -4.,  0.0, Ly +  ly, 0.0,
                     -5.,  0.0, Ly + 0.0,  Lz,
                     -6.,   Lx, Ly + 0.0,  Lz, 
                     -7.,   Lx, Ly +  ly,  Lz, 
                     -8.,  0.0, Ly +  ly,  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[2].Set (/*NDim*/3, /*Tag*/-3, /*NVert*/8,
                     -1.,  0.0, Ly + ly + 0.0, 0.0,
                     -2.,   Lx, Ly + ly + 0.0, 0.0, 
                     -3.,   Lx, Ly + ly +  Ly, 0.0, 
                     -4.,  0.0, Ly + ly +  Ly, 0.0,
                     -5.,  0.0, Ly + ly + 0.0,  Lz,
                     -6.,   Lx, Ly + ly + 0.0,  Lz, 
                     -7.,   Lx, Ly + ly +  Ly,  Lz, 
                     -8.,  0.0, Ly + ly +  Ly,  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);

        // middle
        blks[3].Set (/*NDim*/3, /*Tag*/-4, /*NVert*/8,
                     -1.,  0.0, 0.0, Lz + 0.0,
                     -2.,   Lx, 0.0, Lz + 0.0, 
                     -3.,   Lx,  Ly, Lz + 0.0, 
                     -4.,  0.0,  Ly, Lz + 0.0,
                     -5.,  0.0, 0.0, Lz +  lz,
                     -6.,   Lx, 0.0, Lz +  lz, 
                     -7.,   Lx,  Ly, Lz +  lz, 
                     -8.,  0.0,  Ly, Lz +  lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[4].Set (/*NDim*/3, /*Tag*/-5, /*NVert*/8,
                     -1.,  0.0, Ly + 0.0, Lz + 0.0,
                     -2.,   Lx, Ly + 0.0, Lz + 0.0, 
                     -3.,   Lx, Ly +  ly, Lz + 0.0, 
                     -4.,  0.0, Ly +  ly, Lz + 0.0,
                     -5.,  0.0, Ly + 0.0, Lz +  lz,
                     -6.,   Lx, Ly + 0.0, Lz +  lz, 
                     -7.,   Lx, Ly +  ly, Lz +  lz, 
                     -8.,  0.0, Ly +  ly, Lz +  lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[5].Set (/*NDim*/3, /*Tag*/-6, /*NVert*/8,
                     -1.,  0.0, Ly + ly + 0.0, Lz + 0.0,
                     -2.,   Lx, Ly + ly + 0.0, Lz + 0.0, 
                     -3.,   Lx, Ly + ly +  Ly, Lz + 0.0, 
                     -4.,  0.0, Ly + ly +  Ly, Lz + 0.0,
                     -5.,  0.0, Ly + ly + 0.0, Lz +  lz,
                     -6.,   Lx, Ly + ly + 0.0, Lz +  lz, 
                     -7.,   Lx, Ly + ly +  Ly, Lz +  lz, 
                     -8.,  0.0, Ly + ly +  Ly, Lz +  lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);

        // upper
        blks[6].Set (/*NDim*/3, /*Tag*/-7, /*NVert*/8,
                     -1.,  0.0, 0.0, Lz + lz + 0.0,
                     -2.,   Lx, 0.0, Lz + lz + 0.0, 
                     -3.,   Lx,  Ly, Lz + lz + 0.0, 
                     -4.,  0.0,  Ly, Lz + lz + 0.0,
                     -5.,  0.0, 0.0, Lz + lz +  Lz,
                     -6.,   Lx, 0.0, Lz + lz +  Lz, 
                     -7.,   Lx,  Ly, Lz + lz +  Lz, 
                     -8.,  0.0,  Ly, Lz + lz +  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[7].Set (/*NDim*/3, /*Tag*/-8, /*NVert*/8,
                     -1.,  0.0, Ly + 0.0, Lz + lz + 0.0,
                     -2.,   Lx, Ly + 0.0, Lz + lz + 0.0, 
                     -3.,   Lx, Ly +  ly, Lz + lz + 0.0, 
                     -4.,  0.0, Ly +  ly, Lz + lz + 0.0,
                     -5.,  0.0, Ly + 0.0, Lz + lz +  Lz,
                     -6.,   Lx, Ly + 0.0, Lz + lz +  Lz, 
                     -7.,   Lx, Ly +  ly, Lz + lz +  Lz, 
                     -8.,  0.0, Ly +  ly, Lz + lz +  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);
        blks[8].Set (/*NDim*/3, /*Tag*/-9, /*NVert*/8,
                     -1.,  0.0, Ly + ly + 0.0, Lz + lz + 0.0,
                     -2.,   Lx, Ly + ly + 0.0, Lz + lz + 0.0, 
                     -3.,   Lx, Ly + ly +  Ly, Lz + lz + 0.0, 
                     -4.,  0.0, Ly + ly +  Ly, Lz + lz + 0.0,
                     -5.,  0.0, Ly + ly + 0.0, Lz + lz +  Lz,
                     -6.,   Lx, Ly + ly + 0.0, Lz + lz +  Lz, 
                     -7.,   Lx, Ly + ly +  Ly, Lz + lz +  Lz, 
                     -8.,  0.0, Ly + ly +  Ly, Lz + lz +  Lz,
                     -10.,-20.,-30.,-40.,-50.,-60.);

        size_t Nx=10, Ny=6, Nz=6;
        bool   nonlin = true;
        bool   T      = true;
        double a      = 0.5;
        if (nonlin)
        {
            blks[0].SetNx(Nx);  blks[0].SetNy(Ny, -a, T);  blks[0].SetNz(Nz, -a, T);
            blks[1].SetNx(Nx);  blks[1].SetNy( 1);         blks[1].SetNz(Nz, -a, T);
            blks[2].SetNx(Nx);  blks[2].SetNy(Ny,  a, T);  blks[2].SetNz(Nz, -a, T);
            blks[3].SetNx(Nx);  blks[3].SetNy(Ny, -a, T);  blks[3].SetNz( 1);
            blks[4].SetNx(Nx);  blks[4].SetNy( 1);         blks[4].SetNz( 1);
            blks[5].SetNx(Nx);  blks[5].SetNy(Ny,  a, T);  blks[5].SetNz( 1);
            blks[6].SetNx(Nx);  blks[6].SetNy(Ny, -a, T);  blks[6].SetNz(Nz,  a, T);
            blks[7].SetNx(Nx);  blks[7].SetNy( 1);         blks[7].SetNz(Nz,  a, T);
            blks[8].SetNx(Nx);  blks[8].SetNy(Ny,  a, T);  blks[8].SetNz(Nz,  a, T);
        }
        else
        {
            blks[0].SetNx(Nx);  blks[0].SetNy(Ny);  blks[0].SetNz(Nz);
            blks[1].SetNx(Nx);  blks[1].SetNy( 1);  blks[1].SetNz(Nz);
            blks[2].SetNx(Nx);  blks[2].SetNy(Ny);  blks[2].SetNz(Nz);
            blks[3].SetNx(Nx);  blks[3].SetNy(Ny);  blks[3].SetNz( 1);
            blks[4].SetNx(Nx);  blks[4].SetNy( 1);  blks[4].SetNz( 1);
            blks[5].SetNx(Nx);  blks[5].SetNy(Ny);  blks[5].SetNz( 1);
            blks[6].SetNx(Nx);  blks[6].SetNy(Ny);  blks[6].SetNz(Nz);
            blks[7].SetNx(Nx);  blks[7].SetNy( 1);  blks[7].SetNz(Nz);
            blks[8].SetNx(Nx);  blks[8].SetNy(Ny);  blks[8].SetNz(Nz);
        }

        Mesh::Structured mesh(/*NDim*/3);
        mesh.Generate (blks, /*o2*/false);

        DEM::Domain dom;
        dom.GenFromMesh(mesh,0.1,3.0,true,false);
        dom.WriteXDMF("pullout");

    }

    return 0;
}
MECHSYS_CATCH
